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ABSTRACT 

Eulerian hydrodynamical simulations are a powerful and popular tool for modeling 
fluids in astrophysical systems. In this work, we critically examine recent claims that 
these methods violate Galilean invariance of the Euler equations. We demonstrate 
that Eulerian hydrodynamics methods do converge to a Galilean-invariant solution, 
provided a well-defined convergent solution exists. Specifically, we show that numer- 
ical diffusion, resulting from diffusion-like terms in the discretized hydrodynamical 
equations solved by Eulerian methods, accounts for the effects previously identified as 
evidence for the Galilean non-invariance of these methods. These velocity-dependent 
diffusive terms lead to different results for different bulk velocities when the spatial 
resolution of the simulation is kept fixed, but their effect becomes negligible as the res- 
olution of the simulation is increased to obtain a converged solution. In particular, we 
find that Kelvin-Helmholtz instabilities develop properly in realistic Eulerian calcula- 
tions regardless of the bulk velocity provided the problem is simulated with sufficient 
resolution (a factor of 2-4 increase compared to the case without bulk flows for realistic 
velocities). Our results reiterate that high-resolution Eulerian methods can perform 
well and obtain a convergent solution, even in the presence of highly supersonic bulk 
flows. 

Key words: hydrodynamics-instabilities- mcthods:numerical 



1 INTRODUCTION 

Eulerian methods have been the tool of choice in compu- 
tational fluid dynamics for over five decades. Many suc- 
cess ful Eulerian met hods in popular use descended from 
the iGodunovl dl95Sh scheme that combines the analyti- 
cal Riemann solution of the Euler equation^ with the 
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1 There are also Eulerian astrophysical hydrodynamics codes 
that do not use a God un ov scheme, such as the ZEUS code 
jStone fc Normanlll992alrb1: IClarkelll996l ; lHaves et aljliool) and 
the code bv lRvu" et**al*Tll993l) b ased on the tot al variation dimin- 
ishing flux-corrected method bv lHartenl Jl983T> 



upwind scheme of ICourant et alj (| 19521 ) to numerically 
evolve fluid systems on a discretized mesh. These Godunov- 
type schemes, as such methods are commonly called, have 
been further engineered to include higher-order spatial 
reconstructions of the fluid distr ibution based on piece- 
wise linear (e.g.. Ivan Leerl 1 19771) . parabolic (e.g., PPM, 
IColella fc Woodward!! 19841 ). or, more generally, higher-order 
weighted essent ially non-oscillatory interpolation schemes 
|Liu et al.| [T994) . Eulerian methods have also become quite 
popul ar for addressing problems in Newtonian astrophysics 
(e.g.jFrvxell et all! 19891; ICen et al.lll990l;lBrvan et al.lll994l ; 
Quilis et al.l | l996t[Yepes et al.lll997l : IWada fc Normanlll999l ; 



tive Mesh Refinement (AMR, e.f 


Brvan & Norman! 


1997|; 


Khokhlovl Il998l; iTruelove et all 


19981; iFrvxell et all 


2000|; 


Plewa & Miiller 200ll; iKravtsov et al.ll2002l; Tevssierl 


2002; 


Quilis 2004; Wang et al. 2008). Given their wide-spread use 
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in computational astrophysics, an understanding of the fun- 
damental limitations of such codes is important for inter- 
preting the astrophysics of hydrodynamical systems that 
cannot be accessed through laboratory experiments. 

While Eulerian astrophysical simulation codes rou- 
tinely demonstrate excellent performance on idealized test 
cases, some shortc omings of these methods are known (e.g., 
iQuirkl I1994L 120051 ). Recently, several studies have focused 
on the differences produced by Eulerian codes in refer- 
ence frames moving with different veloc it ies w ith respect to 
the computational grid. IWadslev et al.l (|200&t ) emphasized 
the role of diffusion in altering the developm ent of Kelvin- 
Helm holtz instabilities in the FLASH code (iFrvxell et al 



l200d\ simulations of bouyant, hot bubbles. iTasker et al 



(2008) simulated the advection of otherwise static, self- 
gravitating gas clouds, and showed that the performance 
of FLASH and the PPM version of Enzo (IBrvan fc Normanl 



(IBi 

1997l:lBrvanll999l:lNorman fc Brvanll999l : lBrvan et al.ll200ll : 



Q'Shea et alj 20041 ) in maintaining the centroid and density 
profile of the gas cloud depended on its velocity wit h respect 
to the static computational grid. Most recently, ISpringell 
(2009) motivated the development of the new Lagrangian- 
Eulerian moving-mesh code AREPO by demonstrating that 
with fixed grid Godunov solvers Kelvin-Helmholtz instabili- 
ties may not develop and evolve properly when the interface 
between the two fluids has a large bulk velocity with respect 
to the grid. These apparent failures of Eulerian codes have 
been discussed in terms of "Galilean non-invariance," which 
in this context means that for initial conditions that move 
with different uniform bulk velocities with respect to the 
computational grid but are otherwise identical, numerical 
solutions obtained with Eulerian codes may depend on the 
chosen bulk velocity. 

Given the ubiquity of supersonic bulk motions in as- 
trophysical scenarios, these results are potentially damning 
for the application of stationary mesh Eulerian codes to 
galaxy and structure formation. The purpose of this work 
is to critically examine the performance of Eulerian hy- 
drodynamical codes for simulating systems with supersonic 
bulk motions, and to clarify both the nature and meaning 
of the velocity-dependent differences highlighted in previ- 
ous st udies. Specifically, we use the Eul erian mesh codes 
ART (iKravtsov et all I2002T) and Enzo (IBrvan fc Normanl 
1997l;lBrvanll999l;lNorman fc Brvanll999l ; lBrvan et al.ll200ll ; 
Q'Shea et al.ll2004r ) to simulate the development of Kelvin- 
Helmholtz i nstabilities in t est calculations similar to those 
presented in ISpringell (2009 V We employ statistical measures 
to quantify convergence and error of the calculations in ad- 
dition to an extensive visual comparison of the solutions. We 
show that the effects discussed above are not a consequence 
of Galilean non-invariance of Riemann solvers, but rather 
a result of diffusive errors accumulated during advection of 
fluid through the computational grid. The effects of these 
errors are thus particularly acute in systems where pertur- 
bations and the interface between fluids are under-resolved. 
We demonstrate that with a proper initial setup the Eule- 
rian methods produce a convergent solution at large bulk 
velocity as the resolution of the simulation is increased. 

The paper is organized as follows. In SJH we discuss the 
origin of numerical diffusion in the Eulerian method and 
illustrate its effects using simulations of contact discontinu- 
ities. Readers familiar with the effects of numerical diffu- 



sion should proceed to f|3l where we review previous sim- 
ulations of Kelvin-Helmholtz instabilities and the related 
claims of Galilean non-invariance in Eulerian methods. In 
33] we present a new, better-behaved test calculation of 
Kelvin-Helmholtz instability and study the development of 
the instability over a range of resolutions and supersonic 
bulk motions. We study the statistical and error properties 
of the Kelvin-Helmholtz simulations and use these statistics 
to critically examine the apparent Galilean non-invariance 
of Eulerian simulation codes. We discuss our results in ^5] 
and present a summary in ^6] 



2 NUMERICAL DIFFUSION 

Computational Eulerian hydrodynamical codes calculate the 
evolution of fluid systems using a discretized approximation 
to Euler's equations. When modeling the conservative form 
of Euler's equations, the change in quantities like density 
or energy integrated over cell units of size Ax in the dis- 
cretized mesh over a time step At will correspond to the 
flux of those quantities across the cell boundaries over the 
same time interval. Fluid interactions between cells then fun- 
damentally involve calculations of the fluxes, which can be 
approximated u sing solutions to the Riemann problem (i.e., 
lGodunovlll959r ) o r through other mea ns ( e.g., the flux cor- 
rected methods of Boris fc Booklll973l . and lHarten|[l983l . sec 
also Chapter 21 of lLanevlll998h . Since these numerical ap- 
proximations to the physical fluxes exchanged between fluid 
volumes during the time interval At are discretized, there 
is a truncation error associated with the numerical approxi- 
mation. Missing or extraneous higher order terms in the dis- 
cretized numerical approximation can appear as an effective 
viscosity or thermal conductivity and lead to the smearing 
or dispersion of features in fluid flow. We will refer to smear- 
ing effects as numerical diffusion, while effects that change 
the wave speed of features in the fluid will be labeled nu- 
merical dispersion 0. Clear discussion s about the effects o f 
numerical diffusion can be found in iBoris fc Bookl (I1973D 
and lLanevI i| 19981 ). 

The strength of numerical diffusion depends on the 
method chosen to model fluid systems. Lagrangian methods 
integrate the convective derivative form of the mass con- 
servation equation directly, and therefore suffer from small 
diffusive truncation errors. Eulerian methods calculate the 
advective term in the mass conservation equation explicitly, 
which can lead to an appreciable diffusive truncation error 
upon discretization. Some Eulerian m ethods, such as Flux - 
Corrected Transport algorithms (e.g., IBoris fc Booldll973h . 
include an explicit numerical diffusion term proportional to 
a second spatial derivative that owes to their forced con- 
servative and non-negative properties (the "flux correction" 
refers to the explicit artificial anti-diffusion used to correct 
this truncation error term). 

For Godunov-type methods based on Riemann solvers, 
differences in the amount of numerical diffusion can arise 



2 We note that the numerical diffusion owing to truncation error 
in the Eulerian method is very distinct from artificial viscosity 
employed in Smoothed Particle Hydrodynamics to improve shock 
capturing, and the two should not be confused. 



Computational Eulerian Hydrodynamics and Galilean Invariance 3 



from the approximations made in constructing the dis- 
crete representation of the local fluid flow on the compu- 
tational mesh. In Godunov-type methods, the numerical 
flux between cells is determined by the known solution of 
the piecewise-constant Riemann problem. The resulting flux 
across the cell face is then determined by the properties of 
fluid states on either side of the face, the cell size, and the 
time step size. Resolution determines the region used to av- 
erage the fluid properties for finding the initial states in 
the Riemann problem. The averaging procedure introduces 
numerical diffusion, and can be counter-acted by higher spa- 
tial resolution. Improving the quality of the approximation 
to the fluid states used in the Riemann problem can also 
decrease the amount of numerical diffusion, so the method 
used to model the shape of the fluid flow on the grid can 
change the diffusivity of the method. F or instance, the local 
flow can be approx imated by c onstant (lGodunovlll959f), lin 



apr. 

ear (|van Leerll 19771) , parabolic llColella fc Woodward!! 19841 ). 



or higher-order piecewise polynomial (|Liu et al.lll994l ) inter- 
polations on the discrete mesh. Higher-order interpolations 
improve the local approximations used in reconstructing the 
fluid flow and calculating the initial states to the Riemann 
problem, and therefore will suffer from less numerical dif- 
fusion. In general, the strength of numerical diffusion will 
also depend on the local flow velocity. This velocity depen- 
dence arises because, in the presence of a large advective 
flow, more time steps are used and more local averages are 
performed. Additionally, with a large bulk velocity a larger 
(Lagrangian) region of the fluid is averaged to calculate the 
input states for the Riemann problem. 

We can illustrate how numerical diffusion affects the 
shape of the local fluid distribution by simulating the advec- 
tion of contact discontinuities. In the absence of numerical 
diffusion, the square wave should be perfectly advected and 
the contact discontinuities would remain sharp. However, 
as these simple tests will illustrate, numerical diffusion will 
act to soften the contact discontinuities in a resolution- and 
velocity-dependent manner. The effects of numerical diffu- 
sion in these tests will prove to be informative for simula- 
tions of the Kelvin-Helmholtz instability in Sj3]and[4] 

For the simple problem of the advection of a waveform 
with a constant velocity v, the advected quantity p (e.g., the 
density) obeys the partial differential advection equation 



dp ^dp _ q 
dt dx 



(1) 



as a function of position x and time t. The solution of this 
equation is simply p(x, t) = p(x — vt, 0), as the initial wave- 
form advects with a constant ve locity v. 

However, as discussed by iTorol (|l997l . see, e.g., his 
§5.2.1), numerical methods for solving the advection equa- 
tion (or Euler's equations) actually solve a slightly mod- 
ified equation (to some approximate order). For example, 
in the case o f the simple first-order upwind scheme of 
ICourant et al.l (|l952l ) the modified equation solved by the 
numerical method is 



(2) 



dt dx dx 2 
In this advection-diffusion equation the right-hand side acts 
as a form of numerical diffusion with a diffusion constant a. 
The solution of Equation [2] will differ from the solution of 
Equation[T]if and will be characterized by a progres- 



sive smearing of the original waveform. The detailed behav- 
ior of the solution to Equation [2] will then depend on the 

diffusion constant a. 

F or the first-order upwind scheme of ICourant et al.l 
(1952), the diffusion constant is 



a = -vAx(l - |c|), 



(3) 



where c = vAt/ Ax is the Courant number (numerical sta- 
bility requires |c| ^ 1), At is the timestep, and Ax is the 
spatial grid size. One then expects that the diffusive error 
induced through the truncation of Equation [1] into Equa- 
tion [2] by the discretization of the numerical scheme will 
decrease with the grid size but increase with the advective 
velocity. As Ax — > or v — > 0, the pure advection equa- 
tion Jl} is recovered. Higher-order Eulerian methods (such 
as those used in this paper) can change the form of Equa- 
tions [2] or [3] and also introduce dispersive terms that scale 
as high-order odd-power spatial derivatives. However, as we 
will show, in higher-order methods the strength of numerical 
diffusion will still increase with advection velocity and de- 
crease with increasing spatial resolution. In the remainder of 
this section, we will use square wave advection simulations 
to illustrate these numerical features of Eulerian methods. 

Unless otherwise noted, the simulations presented in 
this paper use the Eulerian code ART with piecewise- linear 
recon struction and an exact Riemann solver (IColella fc Glad 
Il985[). based on the adaptive refinement strategy developed 
bv lKhokhlovl (| 1998h . For the following square wave advection 
simulations, time steps were determined using the Courant- 
Friedrichs-Lewy condition with a parameter cfl = 0.6. ART 
uses a dual-energy formulation similar to that of lBrvan et al.l 
(1995), such that the internal energy equation is followed 
separately when the local flow is kinetic-energy dominated 
and effectively pressureless. We have checked that similar 
results can be obtained using an entropy equation instead 
of the int ernal energy in the dua l-energy f o rmula tion, as dis- 
cussed bv lRvu et ail (|l993h and lSpringej (|2009l ). 

The one-dimensional square wave density is initialized 
with p = 5 for positions \x — 0.5| ^ 0.25 and p = 1 for 
\x — 0.5 1 > 0.25. The system has a constant pressure P — 1 
and an adiabatic index 7 = 5/3. In a first set of tests, the 
wave is advected to the right with a velocity v = 10 in a 
periodic box such that the wave travels through the box 
ten times over the final simulation time t — 1. To illustrate 
the role of resolution on the strength of numerical diffu- 
sion, the simulation is performed with grid resolutions of 
N = [64,128,256,512]. The left panel of Figure [T] shows 
the final square wave density distribution at time t — 1 as a 
function of resolution, compared with the initial distribution 
(thin solid line) . At low resolution (N = 64, red dashed dot- 
ted line), numerical diffusion smears out each contact discon- 
tinuity over approximately twelve cells, or roughly ~ 20% of 
the computational volume. As the resolution increases, the 
contact diffuses out over more cells (~ 18 cells for N = 512, 
black solid line) but less of the computational volume (~ 4% 
for N = 512). The contact is physically better resolved with 
increasing grid size and the diffusive error reduced. The right 
panel illustrates the additional error induced by increasing 
the velocity by a factor of 100 (for N = 512, blue dotted 
line). The highest-resolution simulation has an increased er- 
ror that degrades the effective resolution of the simulation 
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by a factor ~ 4 (comparable to the TV = 128 simulation with 
v = 10 shown in the left panel). 

In these one dimensional simulations, the diffusive error 
can be mitigated through the use of an artificial compression 
(AC) technique. Similar in spirit to explicit anti-diffusion 
terms added to the flux corrected methods, AC simply in- 
creases the slopes used in the reconstruction of the fluid 
on the mesh near contact discontinuities. This approach re- 
duces the effective second order truncation error by limiting 
the influence of the outer cells in the co mputational stencil. 
We use the slo pe steep e ning a pproach of lYand l|l990T ). as im- 
plemented by [Balsaral (11998|) for l inear reconstruction, and 
refer the reader to §2.2 of Balsaral (1 19981 ) for details (see also 



iFrvxell et ail [2000. for an implementation of slope steepen- 
ers for PPM). The right panel of Figure [T] shows the results 
for a TV = 64 grid with bulk velocity v = 10 and AC (solid 
black line). With AC the TV = 64, v = 10 results improve 
to be comparable to the TV = 128, v = 10 results without 
AC. If the same N = 64 simulation with AC is performed 
but with the bulk velocity increased to v — 10 4 (or Mach 
M « 15,000, green dashed line), the results are striking. 
Remarkably, with artificial compression the diffusive error 
becomes almost independent of the bulk velocity and the 
TV = 64, v = 10 4 simulation has almost the same diffusive 
error as the N = 64, v = 10 simulation (and is superior to 
the TV = 512, v — 100 simulation at lOx smaller bulk veloc- 
ity), even as the -/V = 64, v = 10 4 simulation has traversed 
the computational volume 10 4 times and executed 1.2 x 10 6 
individual timesteps. However, a dispersive error has been 
introduced that changes the square wave period by 8% and 
that we have removed in Figure [T] While this fractional dis- 
persive error is only 0.08/10 4 ~ 10 -5 , the error grows to an 
appreciable fraction of the period by t = 1 (as discussed by 
iBoris fc Book|[l973l . this dispersive error may also depend 
on the frequency of the wave form). 

These simulations demonstrate the salient effects of nu- 
merical diffusion on the properties of fluid distributions 
simulated with Eulerian codes. Diffusion limits the sharp- 
ness of fluid distributions, and the averaging of fluid prop- 
erties within cells does not preserve local discontinuities. 
The effects of diffusion can be mitigated through the use 
of higher spatial resolution to improve the local reconstruc- 
tion of the fluid distribution, or through intrinsically less- 
diffusive methods. The presence of a bulk advective velocity 
in the fluid also increases numerical diffusion by increasing 
the number of time steps and local averages of the fluid 
distribution, and can degrade the effective resolution of the 
computational grid. However, the simulation results natu- 
rally improve with increasing grid resolution. With these 
diffusive properties of Eulerian simulations in mind, we will 
now examine simulations of the development of fluid insta- 
bilities in shearing flows. 



3 THE KELVIN-HELMHOLTZ INSTABILITY 

3.1 Kelvin-Helmholtz Instability with a Sharp 
Interface 
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The Kelvin- Helmholtz (KH) instability jHelmholtj Il868l; 
Kelvin! [l910l . see especially Chapter XI of IChandrasekhan 
196lJ) is the unstable growth of perturbations at the inter- 
face between two fluid flows driven by shearing motions. 
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Figure 1. Simulations of a square wave contact discontinuity 
advected with a constant velocity. The contact discontinuity is 
initialized with a density p = 4 for positions \x — 0.5| ^ 0.25 
and p = 1 for \x — 0.5| > 0.25, and a constant pressure P = 1 
(thin black line). The left panel shows the contact discontinuity 
advected with velocity v = 10 simulated with resolutions TV = 64 
(red dash-dotted line), TV = 128 (green dashed line), n = 256 
(blue dotted line), and TV = 512 (thick black line) after time 
t = 1. The numerical diffusive error increases with decreasing 
resolution, and tends to smear out the contact discontinuities. 
The right panel shows the same square wave advected for a time 
t = 1 using a resolution of TV = 64 with advective velocities of 
v — 10 (thick black line) and v = 10 4 (dashed green line), but 
including a r tificial compr ession in the form of slope steepeners 
( Yang 1990; lBalsaralll998h . A phase error of 8.5% in the TV = 64, 
v = 10 4 has been corrected. Artificial compression limits makes 
numerical diffusive error roughly independent of velocity, even for 
advective flows with Mach number M f» 15, 000. Also shown is 
the v = 10 4 simulation with TV = 512 (blue dotted line), which 
has been completely smeared away by diffusive error. 



Perturbations between these fluid phases that grow and be- 
come unstable typically form waves that crest owing to the 
shearing motion in the fluid. The kinetic energy of the shear- 
ing motion powers the instability, and larger shear velocity 
gradients typically increase the proclivity for instabilities to 
develop. In the absence of viscosity and gravity, only inertia 
can exert a stabilizing influence on perturbations and damp 
oscillations before growth commences. 

Numerical simulations of the KH instability previously 
studied in astrophysical contexts incl ude the stability of 
interstellar cloud s in a shearing flow (jMurrav et al.l 1 19931 ; 
IVietri et ai]|l997l : lAgertz et ail 120071 ). t he stripping of gas 
from galaxies by an in tercluster medium (|Quilis et alj|2000l ; 
iMori fc Burkertll2000l ). the formation and ionization state of 
the Magellanic Stream (Bland-Hawthorn ct al. 2007), and 
the s urvivability of high- velocity clouds (|Heitsch fc Putmanl 
2009). Here, we focus on numerical experiments of the KH 
instability in an idealized setting for testing the performance 
of static mesh Eulerian codes in the presence of bulk flows, 
but our conclusions will weigh on the validity of the results 
of many such astrophysical studies. 

A common choice for the initial inhomogeneity that 
gives rise to the KH instability is two uniform fluids sep- 
arated by a surface where the density and shearing veloci- 
ties change discontinuously. The KH instability arising from 
perturbat ions ab o ut th ese i nitial conditions is stu died in 
detail by iKelvinl (|l910h and IChandrasekharl (|l96ll . §100). 
For a surface discontinuity, the growth of any perturbations 
about the surface can be calculated from the Euler equations 
by separating the solution into normal modes. As discussed 
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Figure 2. Kelvin-Helmholtz (KH) instability simulation initial 
conditions for the density (left panel) and ^-direction shear ve- 
locity (right panel) as a function of {/-position. Shown are the 
initial conditions for the Springel (2009) KH simulations ("ICs 
A", black line), as well as a new KH simulation with a smoothly- 
varying density and shear velocity ("ICs B", dashed gray line). 
Both simulations have additional {/-direction velocity perturba- 
tions to seed the instability (see text). The system has a constant 
pressure P = 2.5 and adiabatic index 7 = 5/3. 



bv IChandrasekha^ |l96l|), instability will occur at a per- 
fectly discontinuous interface regardless of the magnitude of 
the shearing velocity. For such initial conditions, this result 
holds generally for some minimum wavenumber, and in the 
absence of gravity or surface tension applies to all wavenum- 
ber s. 

The discontinuous, two-fluid KH instability has been 
used as a test simulation in recent yea rs for evaluat i ng the 
performance of hydrodynamical codes. lAgertz et alJ (|2007l ) 
studied the relative performance of smoothed particle hydro- 
dynamics and Eulerian grid codes in calculating the devel- 
opment of KH instabilities from two-fluid initial conditions. 
ISpringell (2009) also studied KH instabilities in shearing, 
sharp interface between two fluids in two-dimensional sim- 
ulations to test the performance of the Eulerian scheme in 
the presence of bulk flow s. The discontin uous, two-fluid KH 
instability simulations of ISpringell (|2009T ) were performed in 
a unit computational volume in the x — y plane with peri- 
odic boundaries. The initial conditions consisted of a cen- 
tral fluid slab at \y — 0.5 1 < 0.25 with density p\ = 2 and 
v\ = 0.5 (Mach M = 0.35) surrounded by a second fluid at 
\y - 0.5| > 0.25 with density p 2 = 1 and v 2 = -0.5 (Mach 
M = 0.25). The fluids were initialized in pressure equilib- 
rium with P = 2.5 and an adiabatic index of 7 = 5/3. 
A sinusoidal velocity perturbation in the {/-direction of the 
form 



v v {x,y) 



wo sin(n7ra;) x 

(y-0.25) 2 , (y-0.75) 2 



(4) 



cxp 



2<T 2 



+ ■ 



2a 2 



with parameters n = 4, wo — 0.1, and a — 0.05 was added 
to provide a seed for the instability. For reference, Figure 
[2] shows the density and shearing velocity initial conditions. 
We will refer to these discontinuous, two-fluid KH instabil- 
ity initial conditions as "ICs A" . When discussing the Mach 
number of any bulk motions, we will refer to the Mach num- 
ber relative to the sound speed in the dense fluid unless 
oth erwise noted. 

ISpringeil (|2009l ) evolved the system for a time t = 2 
using the new moving-mesh code AREPO using an exact 



Riemann solver i|Torolll997r i in a fixed-mesh mode , and with 
the Eulerian PPM code Athena (IStone et al.ll200Sl ) using the 
linearized solver of iRoe fc Pikei ( 19841 ). The KH simulations 
presented in this paper use the Eulerian code ART, with 
the method described in [J2l un less otherwise noted. As is 
customary, the ART code uses IStrand |l968) dimensional 
splitting to numerically integrate the multidimensional Eu- 
ler equa tions, but we have check ed that using the unsplit 
solver of lGardiner fc Stond (2008) produces similar results. 
For the presented KH simulations, time steps were deter- 
mined using the Courant-Friedrichs-Lewy condition with a 
parameter cfl = 0.6, except near snapshot times where time 
steps were determined by requiring a simulation output ev- 
ery At = 0.01 time interval. 

Figure [3] shows the temporal evolution of this simula- 
tion calculated on a fixed mesh of size 64 x 64 (first row), 
128 x 128 (second row) ,256x256 (third row) , and 2560 x 2560 
(fourth row) at times t = 0.5 (first column), t = 1.0 (second 
column), t = 1.5 (third column), and t = 2.0 (fourth col- 
umn) . A KH instability develops in each simulation, but the 
detailed structure of the growing instability differs between 
simulations with different grid size. The dominant structure 
in the KH instability is the n = 4 mode seeded by the initial 
perturbations (following Equation [4]) . However, a secondary 
set of small-scale eddies that have not been seeded in the ini- 
tial conditions also develop. The development of these small- 
scale instabilities that increase in complexity with increas- 
ing numerical resolution can be directly related to the cho- 
sen fluid interface. As noted bv IChandrasekhail <|l96lh . the 
discontinuous density and shearing velocity distributions of 
the initial conditions allow for perturbations of all wavenum- 
bers to be unstable to growth. An increase in the resolution 
broadens the range of unstable wavelengths available for ex- 
citation by, e.g., secondary waves generated by the seeded 
n — 4 instability or numerical noise, and we further discuss 
these mechanisms below. 

As demonstrated bv ISpringell l|2009h . the simulation of 
the initial conditions ICs A changes dramatically if a uni- 
form bulk flow is added to the fluid. Figure [4] shows the 
results of the simulation at t = 2 of ICs A with a bulk flow 
of v = 10 (Mach M = 6.9) in the {/-direction with resolu- 
tions of N — 64 (left panel), N = 128 (middle panel), and 
N = 256. Each of these panels can be compared directly with 
the results at t = 2 for th e same re s olutio n in Figure [3] and 
are clearly quite different. ISpringeil (|2009l ) states that these 
differences are "direct evidence for a violation of Galilean 
invariance of the Eulerian approach." Although our results 
clearly confirm that the KH instability does not develop in 
the N = 64 simulation, Figure [3] shows that the instability 
does develop at higher resolution and hints at a convergence 
toward a single prominent n = 4 mode instability. 

While qualitative differences between the results of the 
KH instability simulation using initial conditions ICs A are 
apparent in Figures [3] and 2] a quantitative comparison 
would be preferable. A common characterization of a simu- 
lation with a known solution is the error norm, such as the 
L\ error norm given by 



1 



(5) 



where N is the number of computational cells, ft is a prop- 
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Figure 3. Kclvin-Hclmholtz instability simulation of ICs A. Shown is the temporal evolution of the simulation with a mesh resolution 
of 64 X 64 (first row), 128 X 128 (second row), 256 X 256 (third row), and 2560 X 2560 (fourth row) at times t = 0.5 (first column), t = 1.0 
(second column), t = 1.5 (third column), and t = 2.0 (fourth column). 



crty of the ith cell, and /true is the "true" property of same 
cell in the known solution, and the summation runs over all 
N cells. Unfortunately, error norms are useless for evaluat- 
ing the KH simulation of ICs A because there is no con- 
vergence with increasing N and no known solution. How- 



ever, other useful statistical measures can be constructed to 
provide a quantitative gauge of the qualitative differences. 
For instance, the global correlation of the simulations at 
fixed time could be compared using, e.g., Pearson's product- 
moment coefficient. However, the instabilities develop over 
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Figure 4. Kelvin-Helmholtz instability simulation of ICs A including a uniform bulk flow of v = 10 (Mach M = 6.9) in the ^-direction. 
Shown is the computational grid at time t = 2, corresponding to ten full advections of the fluid through the box. The simulation was 
performed on a mesh with TV = 64 (left panel), TV = 128 (middle panel), and TV = 256 (right panel) cells on a side. These results can 
be compared directly with the simulation results at t = 2 shown in Figure [3] Note that the instability fails to develop with resolution 
TV = 64. 



a limited range of the computational volume and simula- 
tions of ICs A with qualitatively very different development 
of the instability would be highly correlated. A more use- 
ful, targeted statistic would track the growth and amount 
of mixing in the instability, but discriminate between the 
sharp features present in the simulations of Figure [3] and 
the diffusive features present in Figure [4] 

After some experimentation, a unit-free measure of the 
root-mean-squared fluctuations in the simulation at fixed y- 
position was found to provide a useful description of insta- 
bility growth and complexity. For a property /, the average 
(/) and variance for each row in the computational mesh 
is calculated. The ratio cr//(/) is then averaged over the 
computational volume as 



W(/> 



dyo f (y)/(f(y)) 



dy 



(6) 



We will refer to the quantity defined by Equation [6] as the 
" mixing statistic" . Ana logous mixing measures were used 
bv lWadslev et ail l|2008l ). Throughout the rest of the paper, 
when the L\ error norm or mixing statistic are used to com- 
pare simulations of differing resolutions, the simulations are 
rebinned to the minimum resolution (usually TV = 64) using 
the IDL function "CONGRID " with the cubic interpolatio n 
value set to CUBIC = -0.5 (|Park fc SchowengerdtJ Il983h . 
Simulations with bulk flows are shifted to align with the 
computational grid as if the measurements were performed 
in the moving frame. 

Figure [5] shows the mixing statistic for the KH insta- 
bility simulation of ICs A as a function of time. Shown are 
the mixing statistics for the density p and the entropy func- 
tion s = P/p 1 for simulations with N = 64, N = 128, 
and TV = 256 both with and without a bulk flow velocity 
of v = 10 (Mach M — 6.9) in the ^-direction. The mixing 
statistic quantifies the qualitative impression that less mix- 
ing occurs in the simulations with a large bulk flow, and 
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Figure 5. "Mixing statistic" for the Kelvin-Helmholtz instability 
simulation of ICs A with time (see Equation [6} . Shown is a di- 
mensionless measure of the root-mean-squared density (left panel) 
and entropy (right panel) fluctuations from the growing KH in- 
stabilities in simulations with TV = 64 (black), TV = 128 (blue), 
and TV = 256 (red) resolution and without bulk flows, as well as 
TV = 64 (black dashed), TV = 128 (blue dashed), and TV = 256 (red 
dashed) simulations with v = 10 (Mach M = 6.9) velocity bulk 
flow along the y-direction. In the simulations without bulk flows, 
the instabilities grow at different rates and by different amounts. 
In the simulations with bulk flows, the instabilities become bet- 
ter defined with increasing resolution but have yet to converge at 
TV = 256 resolution. The instability mostly fails to develop for the 
lowest resolution (TV = 64) simulation with v = 10 bulk flows, as 
noted by Springel (2009). 



that the vertical extent of the instabilities is less than in 
the simulations without a bulk flow. The instabilities grow 
at different rates depending on the resolution, which occurs 
because perturbations with different wavenumbers k grow 
at different rates. The characteristic KH instability growth 
time scales as r oc fc -1 , so the instabilities in the highest res- 
olution simulation (with larger available wavenumbers) grow 
the fastest. In the TV = 128 and TV = 256 simulations with- 
out bulk flows, the primary n = 4 instabilities seeded in the 
initial conditions actually crest and meld with smaller scale 
instabilities resulting from interactions with waves that have 
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traveled across the computational volume. The cresting of 
the waves corresponds to the decrease in the mixing statis- 
tic for the TV = 128 and N = 256 near time t = 1. For the 
simulations with bulk flows, the instability is greatly sup- 
pressed for N = 64 but does develop at higher resolutions. 
In the N = 128 and N = 256 simulations with bulk flows 
the development of the instability is somewhat slower than 
in the simulations without the bulk flows, but faster than in 
the static TV = 64 simulation. 

Given these results, one may ask is there a correct solu- 
tion to which the simulations should converge with increas- 
ing resolution for any bulk velocity? The simulations shown 
in Figure [3] cover a factor of 40 in resolution, and each in- 
crease in resolution is followed by a corresponding increase 
in the complexity of the small-scale structure of the insta- 
bility. Details of the structures, however, are quite different 
in each case as is their overall evolution shown in Fig. [6] As 
such, the solution does not converge with increasing reso- 
lution to any well-defined configuration in the simulations 
with this setup. This result is not surprising as the initial 
conditions with the sharp interface allow all perturbation 
mode s, both real and numerical, to grow |Chandrasekharl 
Il96ll) . The modes excited by wave interactions or seeded by 
numerical noise depend on the actual numerical resolution of 
the simulation and will be different at different resolutions. 
This result is true not only for the static mesh Eulerian cal- 
culations but also f or the cal c ulatio ns with the moving mesh 
code presented by ISpringell (|2009h . We therefore conclude 
that the system ICs A cannot reliably be used for conver- 
gence studies or for the tests studying the development of a 
KH instability in the presence of a uniform bulk flow. 

3.2 Why Does The Simulation Evolution Depend 
On Bulk Velocity? 

The change in the evolution in the presence of a bulk flow has 
been characterized as e vidence for Ga lilean non-invariance of 
the Eulerian methods (|Springelll2009t ). However, the cause of 
the differences has not been unambiguously identified. First, 
as stated above, the simulations without bulk flows do not 
converge with increasing resolution. The cause of this lack 
of convergence is the excitation of small-scale modes by sec- 
ondary waves driven by the initially seeded n — 4 perturba- 
tion. These waves travel through the low-density fluid, cross 
the computational volume, and interact with the dense fluid. 
The interaction between the waves and the dense fluid drives 
high-frequency oscillations that quickly become unstable. At 
high resolution, numerical noise can contribute additional 
small scale structure to these perturbations. These high- 
frequency modes can become unstable owing to the sharp 
transition between the two fluids in the initial conditions. 
If these small-scale instabilities were suppressed, only the 
initially seeded n — 4 mode would grow. 

In the simulations with a bulk flow, the sharp transition 
between the two fluids is smeared owing to diffusive errors 
generated as fluid is advected through the grid — an in- 
herent property of all Eulerian schemes, incl uding those not 
based on Riemann solvers. As discussed by Chandrasckhar 
|l96ll ). the stability of Kelvin-Helmholtz perturbations of 
different wavenumbers depends strongly on the density and 
shearing velocity gradient present between the two fluids. 
While the connection between the instability of a given mode 



and the nature of the gradient can be extremely complicated, 
as a rule of thumb in the absence of gravity and surface ten- 
sion shallower gradients lead to an effective maximum un- 
stable wavenumber of order the inverse of th e spatial scale of 
the gradient (see the discussion in §102 of I Chandrasekharl 
Il96lf ). The numerical diffusion in the Eulerian scheme is 
strong in the simulations with a large bulk flow and simply 
imposes stability on small-scale perturbations. In the simu- 
lations with large bulk flows, only the initially seeded n — 4 
mode grows with time. The lowest resolution (N = 64) sim- 
ulation with a bulk flow has strong enough numerical diffu- 
sion that the n = 4 is not well resolved and diffuses away 
before the shearing flow can cause the wave to crest. 

The results of these Kelvin-Helmholtz simulations sug- 
gest that numerical diffusion leads to change in the available 
modes that can grow into instabilities for the chosen ini- 
tial conditions. The approximation of the physical laws does 
not explicitly change, but the error induced by numerical 
diffusion simply alters the physical system being modeled. 
While this new interpretation of the origin of the differences 
in this Kelvin-Helmholtz simulation is straightforward, it is 
unwieldly to test in this case because of the somewhat patho- 
logical choice of initial conditions. If the advection-related 
diffusion is the origin of the "Galilean non-invariance" of 
the Eulerian schemes, then the error of a numerical solution 
will depend on the bulk velocity (because the integration to 
a given time will be carried out with more time steps), but 
should decrease with increasing resolution. We discuss this 
behavior further in § [5] Since the error norm of ICs A is 
ill-defined in the case without bulk flows a useful error anal- 
ysis would be difficult. We will need to choose a different 
set of initial conditions for a detailed error analysis of KH 
instabilities in the presence of strong bulk flows. 



4 THE KELVIN-HELMHOLTZ INSTABILITY 
WITH A GRADUAL INTERFACE 

As discussed in fj3] the numerical study of Kelvin-Helmholtz 
instabilities can be complicated by the choice of initial con- 
ditions. If large wavenumber perturbations are unstable, the 
growth of the instabilities can be strongly influenced by the 
development of small scale modes seeded or affected by res- 
olution. As a result, the simulation may not converge with 
increasing resolution. A reasonable solution to this problem 
is to choose the initial conditions such that the stratifica- 
tion of the two fluids is not sharp but gradual. Such a setup 
approximates the interfaces that can arise in simulations of 
real astrophysical systems where boundaries between fluids 
are not perfectly discontinuous. 

We therefore alter t he Kelvin-Helm holtz initial condi- 
tions from those used bv lSpringell (120091 ) through the use of 
a "ramp" function 

1 1 

R ^ V ' = 1 + exp[2(y - 0.25)/A y ] 1 + exp[2(0.75 - y)/A v ] ' ^ 

where we set the parameter A y — 0.05. The new density 
distribution is initialized to 

p(y) = pi + R{y)[ P 2 - pi], (8) 

and the shearing velocity distribution is changed to 

v x {y) = vi+R(v)[v2-vi]. (9) 
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Figure 6. Kclvin-Hclmholtz instability simulation of ICs B. Shown is the temporal evolution of the simulation with a mesh resolution 
of 64 X 64 (first row), 128 X 128 (second row), 256 X 256 (third row), and 2560 X 2560 (fourth row) at times t = 0.5 (first column), t = 1.0 
(second column), t = 1.5 (third column), and t = 2.0 (fourth column). 



The parameter values remain p\ = 2, p2 = 1, «i =0.5, and 
t>2 = 0.5, with a constant pressure P = 2.5 and adiabatic 
index 7 = 5/3. The initial velocity perturbation is set to 



v y (x) = wo sin(n7rx), 



(10) 



with wo — 0.1 as before and n = 2. A lower frequency per- 
turbation is chosen to minimize the interaction between in- 
stabilities after they become nonlinear, but we have checked 
that our conclusions are not affected by this choice of per- 
turbation (e.g., using n = 4 leads to similar conclusions, see 
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Figure 7. Kclvin-Hclmholtz instability simulation of initial conditions ICs B at time t = 2. Shown is the simulation density distribution 
for grid resolutions of N = 64 (first column), N = 128 (second column), N = 256 (third column), and N = 512 (fourth column). Each 
grid resolution is simulated with bulk flow velocities of v = [0, 1, 3, 10, 30, 100] (Mach M = [0, 0.7, 2.1, 6.9, 21, 69], top-bottom rows). The 
results of the N = 512, v = run are used to define the L\ error norm for this KH instability simulation. 
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Figure 9. "Mixing statistic" for the Kelvin-Helmholtz instability 
simulation of ICs B with time (see Equation [SJ . Shown is a di- 
mensionless measure of the root-mean-squared density (left panel) 
and entropy (right panel) fluctuations from the growing KH in- 
stabilities in simulations with TV = 64 (black), TV = 128 (blue), 
and TV = 256 (red) resolution and without bulk flows, as well as 
TV = 64 (black dashed), TV = 128 (blue dashed), and TV = 256 
(red dashed) simulations with v = 10 (Mach M = 6.9) velocity 
bulk flow along the y-direction. In the simulations without bulk 
flows, the instabilities grow at nearly the same rate. In the sim- 
ulations with v = 10 bulk flows, the instability growth improves 
with increasing resolution and is comparable to the simulations 
with no bulk flow with a resolution TV = 256 or better. 



Figure 8. Kelvin-Helmholtz instability simulation performed 
with the Piecewise Parabolic Method version of the Enzo code 
(O'Shca et al. 2004). Shown are the simulation results for the den- 
sity at time t = 2 for resolution and ^-direction bulk velocities of 
[N,v] = [64,10] (upper left panel), [N,v] = [128,10] (upper right 
panel), [N,v] = [256,10] (lower left panel), and [N,v] = [512,0] 
(lower right panel). The error convergence rate for PPM recon- 
struction, measured relative to the [TV, v] = [512,0] simulation, 
shows the expected improvement over the linear reconstruction 
results (see Figure as less diffusive methods should perform 
better in the presence of large advective flows. 



fj5]for a discussion). We will refer to these initial conditions 
as "ICs B" , and the corresponding density and shearing ve- 
locity distributions are compared with ICs A in Figure [2] 

The inclusion of a finite gradient in the density and ve- 
locity distribution in ICs B leads to a dramatic suppression 
of small scale features in the growing KH instability. Figure 
[5] shows the temporal evolution of the KH instability arising 
from ICs B with no bulk flow, simulated with grid resolutions 
of TV = 64 (first row), TV = 128 (second row), TV = 256 (third 
row), and TV = 2560 (fourth row). The density distribution 
of the computational volume is plotted at times t = 0.5 (first 
column), t = 1.0 (second column), t = 1.5 (third column), 
and t = 2.0 (fourth column), and is directly comparable to 
the simulations of ICs A shown in Figure [3] The evolution 
of the KH instability is completely dominated by the seeded 
n — 2 perturbation. As with ICs A, the velocity perturba- 
tion drives secondary waves that cross the computational 
volume. These waves travel through the low density fluid 
and collide with the high density fluid after traversing the 
box. In contrast to the evolution of the instability in ICs 
A, these waves do not excite other, higher frequency modes 
in the high density fluid. The transition region between the 
fluids oscillates after interacting with these waves, but the 
density distribution is overstable at large wavenumbers and 
the oscillations damp away. As a result, the evolution of the 
KH instability is nearly independent of the simulation res- 
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Figure 10. L\ error norm for the Kelvin-Helmholtz instabil- 
ity simulation of the initial conditions ICs B with time. Shown 
is the L\ error norm of simulations with no bulk flow and grid 
resolutions TV = 64 (black), TV = 128 (blue), and TV = 256 (red), 
and simulations with bulk flow velocity v = 10 (Mach M = 6.9) 
with resolutions N = 64 (black dashed), TV = 128 (blue dashed), 
N = 256 (red dashed), and TV = 512 (orange dashed). In each 
case, the L\ error norm is measured relative to a TV = 512 sim- 
ulation with no bulk flow. For a bulk flow of velocity v = 10, 
the effective resolution of the simulation is degraded by numeri- 
cal diffusion a factor ~ 4 compared with simulations with no bulk 
flows. 



olution when no bulk flow is included and converges to a 
well-defined solution. 

In an attempt at a comprehensive study of the KH in- 
stability resulting from ICs B, we perform a suite of 24 sim- 
ulations with resolutions TV = [64,128,256,512], with each 
resolution simulation calculated with bulk flow velocities of 
v = [0, 1, 3, 10, 30, 100] (Mach M = [0, 0.7, 2.1, 6.9, 21, 69]) in 
the t/-direction. The simulations were performed in a man- 
ner identical to the simulations of ICs A, with the Courant- 
Freidrichs-Lewy condition parameter cfl — 0.6 and simula- 
tion outputs recorded at time intervals of At = 0.01. 

Figure [7] shows the results of these 24 simulations at 
time t = 2.0, arrayed with resolution increasing to the right 
and bulk flow velocity increasing from v = (top row) to v = 
100 (Mach M = 69, bottom row). The influence of a bulk 
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flow on the evolution of the KH instability growing from ICs 
B is much less dramatic than for ICs A. The diffusive error 
induced by the bulk flow has little influence on the physical 
development of the instability, and only limits the growth 
of the instability for the lowest resolution simulation (TV = 
64) for bulk flow Mach numbers of M > 20 (v > 30). The 
diffusive error clearly decreases and the simulations visually 
appear to converge with increasing resolution at each bulk 
flow velocity. The result is dramatic considering the extreme 
supersonic bulk flow velocities (up to Mach number M ~ 
70) considered. Note that for the run with the largest bulk 
velocity and the highest resolution, the interface has been 
advected through ~ 10 5 computational cells. 

Since the degradation of the computed solution by 
numerical diffusion in the presence of a bulk motion can 
be ameliorated by increasing the resolution, the perfor- 
mance should also improve at fixed resolution when a 
higher-order method is used. The use of PPM reconstruc- 
tion should then result in less diffusion than when lin- 
ear reconstruction is utilized. To test this intuition, we 
use the PPM version of the code Enzo fervan fc Normanl 
1997l;lBrvanll999l;lNorman fc Brvanll999l ; lBrvan et al.ll200ll ; 
Q'Shea et al.ll2004f) . which uses an exact Riemann solver 
and IStrand (|l968T ) dimensional splitting, and perform ex- 
actly the same ramp KH instability simulation. Figure [8] 
shows the results for the density at time t = 2 for sim- 
ulations with resolutions and y-direction bulk velocities of 
N = [64, 128, 256, 512] and v = [10, 10, 10, 0] using the Enzo 
PPM code with cfl = 0.8. First, the results of the Enzo 
PPM and ART simulations are remarkably similar (Figures 
[7] and [8] are directly comparable, with the same color scal- 
ing) . Second, the final density distributions in the lower reso- 
lution (N — 64, 128, 256) simulations with v = 10 bulk flows 
quickly converge with increasing resolution to the reference 
high-resolution ([N,v] — [512,0]) simulation results. As ex- 
pected, at fixed resolution the Enzo PPM results are clearly 
less diffusive than the ART results (e.g., the N = 64, v — 10 
simulation results in Figures [7] and [$} . 

4.1 Error Analysis 

The apparent convergence of the simulation with increasing 
resolution for each bulk flow can be quantified. Statistical 
measures of the instability evolution, including both the L\ 
error norm (Equation [SJ and the mixing statistic (Equa- 
tion O, are well-defined for the ICs B simulation if one sub- 
stitutes the results of a high-resolution simulation for the 
"true" solution. We adopt this approach and define the er- 
ror norm relative to a 512 2 simulation with no bulk flow 
(upper right corner of Figure [Jj), and calculate all statistical 
measures by rebinning simulations to 64 2 resolution when 
necessary. 

For comparison with the results from simulations of ICs 
A, Figure [5] shows the time evolution of the density (left 
panel) and entropy (right panel) mixing statistics for the 
simulation of ICs B with resolutions TV = [64, 128, 256] and 
bulk flow velocities v = [0, 10] (Mach M = [0,6.9]). In dra- 
matic contrast to results of ICs A, the mixing statistic for 
ICs B is roughly independent of resolution for the simula- 
tions with no bulk flow and clearly converges at resolution 
TV = 256 for a bulk flow velocity v = 10. Further, the insta- 
bility grows in both simulations with and without a v = 10 



bulk flow for all the resolutions studied. For simulations with 
a v — 10 bulk flow the growth rate of the instability changes 
with resolution, but the growth rate agrees with the v = 
simulations by resolution N = 256. 

Figure fTOl shows the time evolution of the L\ error norm 
for density (left panel) and entropy (right panel) of the ICs 
B simulations with N = [64, 128, 256] with no bulk velocity 
(v — 0, shown as black, blue, and red lines, respectively) and 
simulations with N = [64,128,256,512] and bulk velocity 
v — 10 (shown as black, blue, red, and orange dashed lines, 
respectively). As expected, the error declines rapidly with 
increasing resolution and increases with increasing time of 
the simulation. For the v = 10 (Mach M = 6.9) bulk flow 
simulations the numerical diffusion degrades the effective 
resolution of the simulation, with the N = 256, v = 10 
simulation performing comparably to the N = 64, v — 
simulation and the N = 512, v = 10 simulation performing 
comparably to the TV = 128, v — simulation. 

The time dependence of the L\ error norm in Figure [TD] 
appears to be self-similar for fixed bulk flow velocity. Fur- 
ther, the logarithmic separation of the L\ error norm for 
simulations with differing resolution appears to be approxi- 
mately constant. With some experimentation, we find that 
the L\ error norm of these simulations scales approximately 
as 

Lx ex A" 2 (l + v) (h55 (l + t)eW64)-° H2,°-°^ (n) 

Figure [11] shows the L\ error norm of all 24 simulations from 
Figure[7Jnormalized by the scaling given by Eg uation l 1 1 1 (left 
panel; the 512 2 , v — simulation error norm is L\ = by 
definition) . The right panel of Figure [11] shows the normal- 
ized Li error norm curves from the left panel divided by the 
normalized L\ error norm for the TV = 256, v — simulation. 
Figure [11] demonstrates that Equation [11] accounts for the 
dependence of the L\ error norm on simulation resolution 
and velocity. At fixed time, the L\ error norm dependence on 
resolution scales as L\ oc 7V~ 2 as expected for the spatially 
second-order accurate Eulerian method used by ART. The 
velocity dependence of the error norm is Li oc (1 + u) ' 55 , 
which is remarkably similar to the L\ oc u ' 5 expected from 
numerical diffusion. We therefore suggest that the velocity- 
dependence of the Li error norm at fixed time is consistent 
with numerical diffusion alone. The time-dependence of the 
Li error norm scaling has only a very weak apparent de- 
pendence on velocity L\ oc (1 + t) 2v , and a moderate 
dependence on the resolution as L\ oc (1 + t) 2 ' JV/ ' 64 - ) . For a 
given bulk flow velocity, Equation QT] can be used to calcu- 
late the necessary resolution required in simulations of this 
KH instabilitj0 to reach an equivalent error for the same 
simulation without a bulk flow. For the setup of ICs B, al- 
though results of runs with and without bulk velocity have 
different errors (expectedly, given the different number of 
time steps and different amount of advection), calculations 
converge to the same Galilean invariant result as the reso- 
lution is increased. 

The same error analysis can be performed for the Enzo 
results as for the ART results above. We use the N = 512, 
v — Enzo PPM simulation, rebinned to N = 64 resolution, 

3 The scaling could conceivably be different in other simulation 
problems. 
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Figure 11. Common L\ error norm dependence on resolution 
and velocity. Shown is the L\ error norm for all 24 simulations in 
Figure [7] normalized by their common dependence on resolution 
[L\ oc TV -2 ) and velocity \L\ oc (1 + v) ' 55 } (left panel). This 
functional dependence is common to all the simulations (right 
panel), as demonstrated by dividing out the time-dependence of 
the TV = 256, v = simulation described by Equation [Tl] 

to define the reference solution in Equation [S] The L\ error 
norms are then measured for the N = [64, 128, 256], v = 10 
Enzo PPM simulations. We find that at early times the error 
norm improves with resolution faster than quadratically, as 
expected for a PPM code. Compared with the results from 
the ART code, these error convergence rates are consistent 
with the precision increase afforded by the PPM reconstruc- 
tion. We therefore suggest that the comparison of the ART 
and Enzo PPM results demonstrate that both codes con- 
verge correctly with increasing resolution for their formal 
order. We therefore expect that other Eulerian codes will 
also produce Galilean invariant solutions to within a speci- 
fied error tolerance for similar hydrodynamical problems. 



5 DISCUSSION 

The simulations of KH instabilities presented in this paper 
suggest that the apparent Galilean non-invariance in Eule- 
rian hydrodynamical calculations owes to numerical diffu- 
sion associated with advection of the fluid through the com- 
putational grid. The diffusive errors arise from the trunca- 
tion errors associated with discretization of spatial and tem- 
poral derivatives in the computational method. The tradi- 
tional meaning of Galilean invariance is that the formulation 
of a specified physical law in two different inertial frames is 
related by a Galilean transformation. While the Euler equa- 
tions are Galilean-invariant, discretized approximations to 
the Euler equations are not guaranteed to obey the same 
transformational properties. However, with increasing spa- 
tial resolution simulations converge to the Galilean invariant 
numerical solution (at least for problems that have conver- 
gent solutions). 

In KH instability simulation initial conditions with dis- 
continuously sharp boundaries between shearing fluids, all 
modes down to the cell scale are unstable. Secondary waves 
driven by the seeded perturbation and numerical noise can 
then generate small-scale instabilities that affect the long- 
term evolution of the system. For such initial conditions, 
numerical diffusion can both alter the frequency range of 
unstable wavenumbers and smear out small-scale features 
in the flow. In the presence of a large bulk flow, the in- 
creased numerical diffusion can thereby suppress the growth 



of small-scale instabilities through these effects. In such tests 
Eulerian methods can produce qualitatively different results 
compared with simulations of the same initial conditions 
that do not include a bulk flow. We emphasize that these 
differences arise because the development of both real and 
numerical small-scale perturbations in the flow are affected 
by diffusive errors associated with advection through the 
grid. 

For physical systems in which all of the perturbations 
are well-resolved, our results show that calculations quickly 
converge to a well-defined, Galilean invariant solution. This 
result holds even in the presence of highly supersonic bulk 
motions, although in this case spatial resolution needs to 
scale approximately as N oc v to counteract increased 
diffusive errors owing to advection. Supersonic bulk flows 
therefore pose no serious problem for Eulerian simulations 
as convergence studies - the final arbiter of any quantita- 
tively credible numerical calculation - will allow one to test 
for a converged solution regardless of the magnitude of the 
bulk flow for realistic initial conditions. Although the re- 
quired increase of resolution may seem like a large price to 
payQ adaptive mesh refinement makes it considerably easier 
to achieve such resolution increases in the relevant regions 
around the fluid interface. We have calculated some of the 
tests presented in the paper with adaptive mesh refinement 
using refinement conditions based on density and entropy 
gradients and obtained results similar to those achieved with 
a uniform grid of higher resolution. 

We caution that numerical diffusion may limit what 
physical systems can be accurately modeled by Eulerian 
methods with arbitrary flow velocities if the resolution is 
poor. Since numerical diffusion can alter the physical system 
being modeled (for instance, in the presence of unrealistic 
discontinuities), simulation results in the Eulerian method 
can depend irrevocably on the local flow velocity in fixed, 
low-resolution calculations. Our studies suggest that hydro- 
dynamical systems altered by numerical diffusion (like ICs 
A) are Galilean invariant to within a specified error toler- 
ance, but these calculations would model a different system 
in the effective absence of numerical diffusion (since, e.g., 
the frequency range of unstable modes has changed). The 
failure of KH instabilities to develop in low-resolution sim- 
ulations with large bulk velocities can be directly rectified 
with sufficient resolution. We therefore suggest that knowl- 
edge of the bulk velocity, local density and shear velocity 
gradients, or other gauges of the potential impact of numer- 
ical diffusion be incorporated into resolution criteria and, for 
adaptive mesh codes, into refinement conditions. 

When Eulerian methods are used to model the growth 
of instabilities in the presence of a bulk velocity, care should 
be exercised to ensure that the perturbations of interest are 
not made overstable by numerical diffusion at low resolu- 
tion. For instance, if the perturbation frequency n in Equa- 
tion [TU] was greatly increased such that the perturbation 
wavenumber was close to the maximum unstable wavenum- 
ber permitted by the system, then numerical diffusion could, 



4 Note that the use of moving mesh to handle the bulk flows is 
also expensive, as it requires both large memory per cell and a 
fairly sophisticated algorithm for handling the motion of the mesh 
and regridding. 
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in principle, artificially alter the growth of the perturbation 
if the mode was not properly resolved. We have repeated 
simulations of ICs B with n = 4 and find similar conver- 
gence properties to those reported for n — 2. However, the 
typical error norm at fixed resolution and velocity is larger 
for n — 4 than for n = 2. Sensibly, smearing in the solution 
owing to numerical diffusion will affect small scales more so 
than large scale modes. 



5.1 Limiting Numerical Diffusion 

Techniques exist for reducing numerical diffusion in Eule- 
rian methods, and an appropriate application of such tech- 
niques could help mitigate velocity-dependent behavior of 
calculated solutions. We demonstrated that intrinsically less 
diffusive Eulerian methods, such as PPM, show improved re- 
sults over codes using linear reconstru ction. The ART simu- 
lations presented in this paper use the Ivan Leerl (|l977l ) slope 
limiter, but the use of less-diffusive slope lim iters durin g 
reconstruction, such as the Superbee limiter 
improve results for some physical systems. Artificial com- 
pression in the form of slope steepeners might also help for 
some problems. For instance, our tests presented in Sj2]show 
that the one dimensional advection of square wave contact 
discontinuities can be made insensitive to even ultrasonic 
(Mach M ~ 10 4 ) advective velocities if the the slope steep- 
ening method of lYand (|l990T ) is used. However, when applied 
to multidimensional problems, such as the KH inst ability 
simula tions presented here, we have found that the lYand 
( 1990) slope steepener produces oddly angular features in 
the fluid flow when used with a IStrand l|l968h dimension- 
ally split solver (we have not tested s l ope st eepeners with 
unsplit solvers) . Similarly, the iHartenl l|l989T ) subcell reso- 
lution method produces excellent results in one dimension 
but is unwieldly in multiple dimensions. Other approaches 
for reducing numeri cal diffusion in Eulerian methods are dis- 
cussed at length by iLanevl l| 19981 ). 

The strength of numerical diffusion can differ on the up- 
wind and downwind side of an initially symmetrical wave- 
form being advected using Godunov-type Eulerian codes. 
The results of §3] show clear evidence for larger numerical 
diffusion in the high-to-low density transition than for the 
low-to-high density transition, and for an alteration of the 
piecewise reconstruction of the initially symmetrical wave- 
form by numerical diffusion. It should be noted that some 
reconstruction methods intrinsically do not produce sym- 
metrical piecewise approximations to smooth functions, de- 
pending on the number of samples of the waveform. For in- 
stance, average-quadratic interpolations of sine and square 
waves using essentially nonoscil latory recons truction are not 
always symmetrical (see §9.2 of |LaneyJ[l998|) ■ Numerical dif- 
fusion may depend on the flow velocity and direction in such 
cases. 



5.2 Gravity and the Kelvin-Helmholtz Instability 

While our study has focused on idealized hydrodynami- 
cal systems, real astrophysical systems where the Kelvin- 
Helmholtz instability operates will be influenced by gravita- 
tional fields. Mixing in a stratified fluid in the presence of 
gravity must overcome gravitational potential energy with 



kinetic energy from the shearing motion. The stabililty of 
the fluid against Kelvin-Helmholtz modes can therefore be 
characterized by the Richardson number that measures the 
relative strength of buoyancy in the fluid and inertia sup- 
plied by the shearing velocity gradient dv/dz as 



J : 



9 dp/dz 
p (dv/dz)' 



(12) 



where g is the local gravitational acceleration along the z- 
direction and p and dp/ dz are the local densi ty and den- 
sity gradient in the fluid. IChandrasekharl (|l96ll ) shows that 
since the kinetic energy of the shearing motion powers the 
instability, if the shearing kinetic energy is too small to over- 
come the gravitational potential energy of the fluid then the 
KH instability will not occur. In terms of the Richardson 
number, it is straightforward to show that the correspond- 
ing necessary (but not sufficient) condition for stability is 
J > 1/40 

In the context of modeling astrophysical systems, our 
results suggest that numerical diffusion may act to change 
the Richardson number by altering the local density and 
shear velocity gradients. Since the numerical diffusion will 
smooth density and velocity gradients by similar amounts, 
the effect of numerical diffusion will typically be to increase 
the stability of KH modes. Simulations of systems with 
modes that have J ~ 1/4 should therefore be checked to 
gauge the influence of any artificial stabilization from nu- 
merical diffusion on the global evolution of the system. 



6 SUMMARY 

In this paper we have presented hydrodynamical simulations 
of Kelvin-Helmholtz (KH) instabilities to study the behavior 
of the numerical solution in the presence of uniform bulk mo- 
tion. We perform these calculations to verify and evaluate 
previous claims in the literature that Godunov-type Eule- 
rian mesh calculations are inheren tly Galile a n non -invariant. 

The KH instabilit y study of ISpringell (20091 was per- 
formed using the ART jKravtsov et al.lll997l,l2002ii hydrody- 
namical code over a range of numerical reso lutions and bulk 
flow velocities. We confirm the results of ISpringell (|2009T ) 
that for low-resolution simulations KH instabilities may not 
develop in the presence of bulk flows, but show that such 
instabilities do develop w i th su fficient resolution. We also 
explain why the ISpringe] (|2009h KH instability simulation 
results generally depend on whether a bulk flow is included. 
Diffusion in the presence of a bulk flow softens the sharp dis- 
continuity between fluids in these initial conditions, thereby 
changing the frequency range of unstable modes. As a result, 
the small-scale structure of the solution changes significantly 
with resolution and in the presence of large advective veloc- 
ity. However, we emphasize that this velocity dependence 
owes to numerical diffusion associated with advection of the 
fluid, not because the analytical Riemann solution is some- 
how Galilean non-invariant (it is invariant). 



5 In the absence of gravity J = 0, and stability depends on 
whether individual modes can be excited. In the presence of other 
effects, such as a Coriolis for ce, instability may requi re higher 
Richardson numbers (see, e.g., lG6mez &: O strikcr 20051). 
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We test this explanation by simulating another Kelvin- 
Helmholtz instability that grows from initial conditions that 
have a gradual interface b etween fluids, in contrast to the 
sharp interface used in the ISpringell (|2009l ) KH test. In the 
new KH simulation, the stability of the fluid to growing 
small-scale modes allows the seeded n = 2 perturbation 
to develop unimpeded as a KH instability, and the simu- 
lation demonstrates convergence with increasing resolution. 
Furthermore, the development of the KH instability occurs 
in the same manner for all bulk flows examined (includ- 
ing supersonic bulk flows with Mach numbers M ~ 70). 
An analysis of the L\ error norm suggests that numerical 
diffusion accounts for the entire error budget of the sim- 
ulations that include bulk flows. We support this conclu- 
sion by demonstrating that the intrinsi cally less diffusive 
Piecewise Parabolic Method code Enzo H^v an fe Norman! 
1997l;lBrvanll999l;lNorman fc Brvanll999l ; lBrvan et al.ll200ll ; 
O'Shea et al.l I2004T ) exhibits more rapid convergence than 
codes using linear reconstruction, and that the simulation 
results produced by ART and Enzo are consistent. For this 
KH instability, the Galilean non-invariance can therefore be 
entirely accounted for by numerical diffusion and can be ef- 
fectively eliminated with increasing numerical resolution. 

Our results suggest that physical systems where numer- 
ical diffusion does not significantly alter the frequency range 
of unstable modes, Godunov-type Eulerian methods will be 
Galilean invariant to within a specified numerical error and 
that this error will decrease with increasing resolution. We 
have demonstrated this result explicitly in the case of a KH 
instability, but we suspect our conclusions will generalize to 
other hydrodynamical instabilities. Similar conclusions can 
be drawn from test cal culations by ot her authors (e.g., the 
Gresho vortex tests by ISprin gcl 2009). These results show 
that there is no generic problem of using the Eulerian meth- 
ods for modeling complicated astrophysical systems with 
large bulk flows. However, overcoming diffusive errors will 
require more stringent resolution when modeling systems 
with large bulk fluid motions. 

As with most numerical calculations, convergence stud- 
ies are essential for evaluating how numerical diffusion in- 
troduces velocity-dependent error into the presented simu- 
lations. If a full converge n ce stu dy is difficult or impossible, 
such as for the ISpringe] (|2009h KH instability simulation, 
caution needs to be exercised in interpreting results in the 
presence of high-Mach number flows. 
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